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Abstract 

Non equilibrium Green's function methods are regularly used to calculate current and charge 
densities in nanoscale (both molecular and semiconductor) conductors under bias. This method is 
mainly used for ballistic conduction but may be extended to include inelastic scattering. In this tutorial 
paper the NEGF equations for the current and charge density matrix are derived and explained in a 
hopefully clear way. 

1 Introduction 

Non equilibrium Green's function methods are regularly used to calculate current and charge densities 
in nanoscale (both molecular and semiconductor) conductors under bias. An overview of the theory of 
molecular electronics can be found in Ref. PP and for semiconductor nanoscale devices see Ref. 

The aim of this text is to provide some intuitive explanations of one particle Green's functions in a 
compact form together with derivations of the expressions for the current and the density matrix. It is 
not intended as a complete stand-alone tutorial, but rather as a complement to Ref. ^|21E1E]. 

2 Green's functions 

Discrete Schrodinger equation: 

H\n)=E\n) (1) 

We divide the Hamiltonian and wavefunction of the system into contact (-ffi,2, \il'i,2)) and device {Hd, 
\^d)) subspaces: 

( Hi \ / iV^i) \ / \i,i) \ 
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where ti,2 describes the interaction between device and contacts. In general we have N contacts {Hi^...^n) 
connecting (ti^...^7v) the device Hd to the reservoirs. Here we will assume that the contacts are independent, 
i.e., there are no cross terms (r) between the different contacts. 
We define the Green's function^: 

iE-H)GiE)^I (3) 

2.1 Why do we want to calculate the Green's function? 

The Green's function gives the response of a system to a constant perturbation 
equation: 

H\il;) = + \v) 

The response to this perturbation is: 

{E-H)\^) = -\v) ^ (5) 

IV) - -G{E)\v) (6) 

Why do we need the response to this type of perturbation? Well, it turns out that it's usually easier 
(see next section) to calculate the Green's function than solve the whole eigenvalue problem^ and most 
(all for the one-particle system) properties of the system can be calculated from the Green's function. 
E.g., the wavefunction of the contact (IV'2)) can be calculated if we know the wavefunction on the device 
{\iJd))- From third row of Eq. 13 

iJ2|V'2)+T2|Vd) = E\iP2) ^ (7) 

{E-H2)\i:2) = ^ (8) 

\i>2) = 92{E)T2\^Pd) (9) 

where g2 is the Green's function of the isolated contact 2 {{E — i?2)ff2 = I)- 

It is important to note that since we have an infinite system, we obtain two types of solutions for the 
Green's functions'^, the retarded and the advanced** solutions corresponding to outgoing and incoming 
waves in the contacts. 

Notation: We will denote the retarded Green's function with G and the advanced with (and maybe 
G^ and G^ occasionally). Here, CAPITAL G denotes the full Green's function and its sub-matrices Gi, 
Gd, Gid etc. Lowercase is used for the Green's functions of the isolated subsystems, e.g., {E — H2) §2 = I- 

Note that by using the retarded Green's function of the isolated contact (172) in Eq. IHlwe obtain the 
solution corresponding to a outgoing wave in the contact. Using the advanced Green's function (g^) would 
give the solution corresponding to an incoming wave. 

^Others may (and do) use the opposite sign. 
^Especially for infinite systems. 

■^When the energy coincides with energy band of the contacts there are two solutions corresponding to outgoing or 
incoming waves in the contacts. 

^In practice these two solutions are usually obtained by adding an imaginary part to the energy. By taking the limit to 
zero of the imaginary part one of the two solutions is obtained. If the limit — > 0'^ is taken the retarded solution is found, 
— > 0~ gives the advanced. This can be seen from the Fourier transform of the time dependent Green's function. 



\v) in the Schrodinger 
(4) 
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2.2 Self-Energy 

The reason for calculating the Green's function is that it is easier that solving the Schrodinger equation. 
Also, the Green's function of the device {Gd) can be calculated separately without calculating the whole 
Green's function (G). From the definition of the Green's function we obtain: 

E-Hi -n \ / Gi Gid Gi2 \ / / \ 

-4 E-Hd -rt Gdi Gd Gd2 - / (10) 

-T2 E-H2 J \ G21 G2d G2 J V / / 

Selecting the three equations in the second column: 

{E-H,)Gid~nGd - (11) 
~TlGid + {E-Hd)Gd-TlG2d = I (12) 
{E-H2)G2d-T2Gd - (13) 

We can solve Eqs.^Jand^lfor Gid and G2d- 

Gid - giTiGd (14) 

G2d = g2T2Gd (15) 

substitution into Ea.lT^ gives: 

- AgiTiGd + {E- Hd) Gd - 4g2T2Gd - / (16) 



from which Gd is simple to find: 

Gd = {E-Hd-J:i-^2y^ (17) 

where Ei = rlgiTi and S2 = T2g2T2 are the so called self-energies. 

Loosely one can say that the effect of the contacts on the device is to add the self-energies to the 
device Hamiltonian since when we calculate the Green's function on the device we just calculate the 
Green's function for the effective Hamiltonian i?ctfoctivc — Hd -I- Si -I- S2. However, we should keep in mind 
that we can only do this when we calculate the Green's function. The eigen-values and -vectors of this 
effective Hamiltonian are not quantities we can interpret easily. 

For "normal" contacts, the surface Green's functions gi and 52 used to calculate the self-energies are 
usually calculated using the periodicity of the contacts, this method is described in detail in appendix B 
of Ref. [3] and in section 3 of Ref. P]- 

2.3 The spectral function 

Another important use of the Green's function is the spectral function: 

A^i{G-G^) (18) 

which gives the DOS and all solutions to the Schrodinger equation! 

To see this we first note that for any perturbation \v) we get two solutions (|'0''^) and IV'"^)) to the 
perturbed Schrodinger equation: 

iE-Hm^-\v) (19) 
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from the advanced and retarded Green's functions: 



= -G\v) (20) 
1^)^ = -G^l^) (21) 

The difference of these solutions (IV--") - IV'^) = (G - G'<)\v) = ~iA\v)) is a solution to the Schrodinger 
equation: 

(E - H){\^^') - 1^^)) = {E- H){G - G^)\v) = {I- I)\v) = (22) 

which means that \^) = A\v) is a solution to the Schrodinger equation for any vector \v)\ 

To show that the spectral function actually gives all solutions to the Schrodinger equation is a little 
bit more complicated and we need the expansion of the Green's function in the eigenbasis: 

G= I l^^f (23) 

k 

where the 6 is the small imaginary part (see footnote 13 , |fc)'s are all eigenvectors^ to H with the corre- 
sponding eigenvalues Cfe. Expanding the spectral function in the eigenbasis gives: 



E + iS-H E-i5-H 
1 

E + id — ek E — iS — Ck 



k 

= y\k){k\ (26) 

where 6 is our infinitesimal imaginary part of the energy. Letting 6 go to zero gives: 

A = 27rY,HE-ek)\k){k\ (27) 

k 

(here S{E — ek) is the delta function) which can be seen since goes to zero everywhere but at 

E = £k, integrating over E (with a test function) gives the 27r5 {E — ek) factor. Eq. 1271 shows that the 
spectral function gives us all solutions to the Schrodinger equation. 



3 Response to an incoming wave 

In the non-equilibrium case, reservoirs with different chemical potentials will inject electrons and occupy 
the states corresponding to incoming waves in the contacts. Therefore, we want to find the solutions 
corresponding to these incoming waves. 

Consider contact 1 isolated from the other contacts and the device. At a certain energy we have 
solutions corresponding to an incoming wave that is totally reflected at the end of the contact. We will 
denote these solutions with IV'i.ri) where 1 is the contact number and n is a quantum number (we may 

^Normalized! 
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have several modes in the contacts). We can find all these solutions from the spectral function ai of the 
isolated contact (as described above). 

Connecting the contacts to the device we can calculate the wavefunction on the whole system caused by 
the incoming wave in contact 1. To do this we note that a wavefunction should be of the form iV'i.n) + 1 V'^) 
where \4'i,n) is the totally reflected wave and \ip^) is the retarded response of the whole system. Putting 
in the anzats \ipi,n) + into the Schrodinger equation gives: 

Hi + n+ \ 

Hd + tI + 4+ (|V^i,„) + |^«)) = E (|^i,„) + \^^)) (28) 
H2 + T2 ) 

E\ll^l,n)+ \ ( H1+T1+ \ 

Tl\^i,n)+ + H, + tI+tI+ I^^) = £; (|V^i,„) + l^"")) (29) 

/ V H2 + T2 J 

H\tp") = i?|^«)-r||V'i,„) (30) 

(note the slight change in notation) and we see that \ip^) is noting else the response of the whole system 
to a perturbation of — Tj^|i/'i,n), c.f., Eq.|51 

IV-^) = Grt (31) 

It is important to realize that the scattering states generated from Eq. 1311 using all possible incoming 
waves from each contact, form a complete^ ON set of solutions to the full Schrodinger equation[S]. Note 
that we have chosen the retarded response which means that the only part of the wave that is traveling 
towards the device is the incoming wave (part of IV'i.n))- We will make full use of this fact below. 

It will be useful to have the expressions for the device wavefunction \ijjd) and contact wavefunction 
(IV'1,2)). The device part is straightforward: 

m = Gdrll^Pi.n) (32) 

and from Eq. Elor Eg. [TBI 

IV'2) = g2T2\lpd) = g2T2GdTl\lpl,n) (33) 

Note, to calculate the wavefunction in the contact containing the incoming wave (contact 1) we need to 
add the incoming wave, giving a slightly more complicated expression: 

l^i) = (l+.9iTiGrfT^) |^i,„) (34) 

Knowing the wavefunctions corresponding to incoming waves in different contacts enables us to fill up 
the different solutions according to the electron reservoirs filling the contacts. 

4 Charge density matrix 

In the non equilibrium case we are often interested in two quantities: the current and the charge density 
matrix. Lets start with the charge density (which allows us to use a self-consistent scheme to describe 
charging) . 

^Except for localized states in the device region. 
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The charge density matrix is defined as: 

p = ^/(fc,M)|^fc)(V'fe| (35) 

k 

where the sum runs over all states with the occupation number f{Ek, /i) (pure density matrix) (note the 
similarity with the spectral function A, in equilibrium you find the density matrix from A and not as 
described below). In our case, the occupation number is determined by the reservoirs filling the incoming 
waves in the contacts such that: ^ 

is the Fermi-Dirac function with the chemical potential {pi) and temperature (T) of the reservoir respon- 
sible for injecting the electrons into the contacts. 

The wavefunction on the device given by an incoming wave in contact 1 (see Eq. I32II is: 

l^d^k) ^ GdAliPi^k) (37) 

Adding up all states from contact 1 gives: 

oo 

Pd[contactl] = J dEj2fiE,f^i)SiE-Ek)\^a.k){^d,k\ (38) 

£=-00 

= J dE f{E,fii)Y,S{E-Ek)GdTl\^Ji,k){^Ji,k\TiGl (39) 
= J dE fiE,fii)GdTl(y2S{E'Ek)\4'i,k){4'i,k\\TiGi (40) 

oo 

= [Eq.|lII= J dE f{E,^i^)Gdrl^TiGl (41) 

£=-00 

introducing the new quantity Fi = r|aiTi = i — we obtain the simple formula: 

oo 

p[from contact 1] = 7^ / dE f{E,fii)GdTiGl (42) 
27r J 



E= 

The total charge density thus becomes a sum over all contacts: 



00 

^^2(fc^ / dEj:fiE.,.)Gdr.Gl (43) 

E 



6 



5 Probability Current 



Having different chemical potentials in the reservoirs filling the contacts gives rise to a current. In the 
next section we will calculate this current in a similar way as the charge density was calculated. But to 
do this we need an expression for the current from the wavefunction. 

In the continuum case we can calculate the current from the velocity operator. However, for a discrete 
Hamiltonian it is not so clear what the velocity operator is. Therefore, we derive an expression for the 
current from the continuity equation (using two contacts). In steady-state, the probability to find an 
electron on the device (^ l^/'ip where the sum runs over the device subspace) is conserved: 



y: ^mm ^ ^ (^^(,1^) + m^'] (44) 



dt ^ dt ^\ dt ' ' dt 

i i 
i 

= ^{{md + ri+T2\ibd)-{MHd + rl+T^^\ij)) (46) 
= J {[{^i\rM - iUrlli^i)] + [{HrM - {U4\^2)]) (47) 

We interpret the term in the first (square) bracket as the incoming probability current into the device 
from contact 1 and the second bracket from contact 2. Generalizing to an arbitrary contact j gives us the 
electric current (at one energy) as the charge (-e) times the probability current: 



ze 



((^,|r,#rf)-(V'd|rj|^,)) (48) 



where ij is defined as positive for a current from the contacts into the device. We can now put in the 
expressions for the wavefunctions in the same way as for the density matrix. 

6 Electrical Current 

To calculate the total current through the device we only need to put in the wavefunction of the device 
and the contacts {\ipd), IV'i); \'>p2)) from Eas. 1321 nU andl?^ and add all the contributions together. Thus the 
current into the device from a incoming wave of one energy (E) in contact 1 {\ipi,n)) through the coupling 
defined by T2 is: 

= -J {{i^2\T2m - {U4\i^2)) (49) 

= -ji{^l,n\TlGy29lT2GdTl\lbi,n) " iTiG^rjgzTzGdrt j^-l^n)) (50) 

= -^(V'l.nlTiGtrt (.gt _ ^2) r2GdTt|Vi,„)) (51) 
= J(V'i,„|TiGtr2Grfr||^i,„) (52) 
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Adding over the modes n and noting that the levels are filled from the reservoir connected to contact 1 
gives (2 for spin): 



oc 

/.frcni = 2^ y" di;/(£;,Mi)5I^(^-^n)(V'i,™|Ti4r2Gdrti^i,„) 



(53) 



E=-cxi 



2e 

'n 

2e 
¥ 

2e 



E= 



E=-oo 



e 



oc 

f dEf{E,n,)Y,S{E-En){il^i,n\Ti\m){m\GlT2GdTl\iJi,n) (54) 

- — oo ' 

J f{E,f,i)Y,{rn\Glr2GaTl(^^6{E-Er,)\i;i,n){i^i,n\\Ti\rn) (55) 

OO 

/ dEf{E,,,^)J2{m\GlT2GdTl^n\m) (56) 

m 

— oo 

oo 

j d£;/(^,Mi)Tr(Gir2Gdri) (57) 

B=-oo 

To get the total current through the device the current from contact two have to be subtracted away: 

oo 

7=-^ J dE{f{E,fx,)-f{E,i^2))Tr(Glr2GdT^) (58) 

E=-oo 

which is exactly the Landauer formula for the current. 
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